{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# file header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys,os,inspect\n",
    "if hasattr(sys.modules[__name__], '__file__'):\n",
    "    _file_name = __file__\n",
    "else:\n",
    "    _file_name = inspect.getfile(inspect.currentframe())\n",
    "CURRENT_FILE_PATH = os.path.dirname(_file_name)\n",
    "sys.path.append(os.getcwd()+\"/../../neuronVis\")\n",
    "from cx_tool import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# tmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = '../../resource/project_amygdala/'\n",
    "create_dir(data_path)\n",
    "sst_samplelist = [\"221058\", \"221297\", \"230058\", \"242655\", \"242656\"]\n",
    "pkc_samplelist = [\"234064\", \"234065\", \"242657\", \"242658\"]\n",
    "crh_samplelist = [\"234171\", \"234172\"]\n",
    "iondata = IONData.IONData()\n",
    "mor_CEA_neuron_list = []\n",
    "sst_CEA_neuron_list = []\n",
    "for s in sst_samplelist:\n",
    "    neuronlist = iondata.getNeuronListBySampleID(s)\n",
    "    for n in neuronlist:\n",
    "        if n['region'][:3] == 'CEA':\n",
    "            sst_CEA_neuron_list.append(n['sampleid']+n['name'])\n",
    "save_scene(sst_CEA_neuron_list,color_pool[0],data_path+'SST')\n",
    "pkc_CEA_neuron_list = []\n",
    "for s in pkc_samplelist:\n",
    "    neuronlist = iondata.getNeuronListBySampleID(s)\n",
    "    for n in neuronlist:\n",
    "        if n['region'][:3] == 'CEA':\n",
    "            pkc_CEA_neuron_list.append(n['sampleid']+n['name'])\n",
    "save_scene(pkc_CEA_neuron_list,color_pool[1],data_path+'PKC')\n",
    "crh_CEA_neuron_list = []\n",
    "for s in crh_samplelist:\n",
    "    neuronlist = iondata.getNeuronListBySampleID(s)\n",
    "    for n in neuronlist:\n",
    "        if n['region'][:3] == 'CEA':\n",
    "            crh_CEA_neuron_list.append(n['sampleid']+n['name'])\n",
    "save_scene(crh_CEA_neuron_list,color_pool[3],data_path+'CRH')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projection neurons from the whole brain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = '../../resource/project_amygdala/'\n",
    "create_dir(data_path)\n",
    "amygdala_heirarchy = {'OLF':['COA','NLOT'],'CTXsp':['LA','BLA','BMA','PA'],'sAMY':['AAA','BA','CEA','IA','MEA']}\n",
    "amygdala_region = []\n",
    "for i in amygdala_heirarchy:\n",
    "    for j in amygdala_heirarchy[i]:\n",
    "        amygdala_region.append(j)\n",
    "if not os.path.exists(data_path+'neuron_information.npy'):\n",
    "    amygdala_projection_neuron_name = []\n",
    "    for i in amygdala_heirarchy:\n",
    "        for j in amygdala_heirarchy[i]:\n",
    "            for k in iondata.getNeuronListByTerminalRegion(j):\n",
    "                if len(k['sampleid'])==6 and k['sampleid'][0]!='0' and len(k['name'])==7:\n",
    "                    neuron = k['sampleid']+k['name']\n",
    "                    if neuron not in amygdala_projection_neuron_name:\n",
    "                        amygdala_projection_neuron_name.append(neuron)\n",
    "    neuron_information = {}\n",
    "    for i in amygdala_projection_neuron_name:\n",
    "        control_console_output(0)\n",
    "        neuron_property = iondata.getNeuronPropertyByID(i[:-7],i[-7:])\n",
    "        neuron_information[i] = {}\n",
    "        neuron_information[i]['cortex_region'] = neuron_property['somaregion']\n",
    "        neuron_information[i]['axonal_projection_terminal_combined'] = get_axonal_statics_by_regions(i,amygdala_region,'terminalregion')\n",
    "        neuron_information[i]['soma_position'] = neuron_property['somapoint']\n",
    "        control_console_output(1)\n",
    "    np.save(data_path+'neuron_information.npy',neuron_information)\n",
    "else:\n",
    "    neuron_information = np.load(data_path+'neuron_information.npy',allow_pickle=True).item()\n",
    "    amygdala_projection_neuron_name = list(neuron_information.keys())\n",
    "if not os.path.exists(data_path+'neuron_terminals.npy'):\n",
    "    neuron_terminals = {}\n",
    "    for neuron in amygdala_projection_neuron_name:\n",
    "        control_console_output(0)\n",
    "        neuron_tree = get_neuron_tree(neuron)\n",
    "        control_console_output(1)\n",
    "        neuron_terminals[neuron] = {'mirror_mark':neuron_tree.root.z>5700,'terminals':[]}\n",
    "        for terminal in neuron_tree.terminals:\n",
    "            if terminal.type==3:\n",
    "                continue\n",
    "            neuron_terminals[neuron]['terminals'].append(terminal.xyz)\n",
    "    np.save(data_path+'neuron_terminals.npy',neuron_terminals)\n",
    "else:\n",
    "    neuron_terminals = np.load(data_path+'neuron_terminals.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(data_path+'neuron_terminals.npy'):\n",
    "    neuron_terminals = {}\n",
    "    for neuron in amygdala_projection_neuron_name:\n",
    "        control_console_output(0)\n",
    "        neuron_tree = get_neuron_tree(neuron)\n",
    "        control_console_output(1)\n",
    "        neuron_terminals[neuron] = {'mirror_mark':neuron_tree.root.z>5700,'terminals':[]}\n",
    "        for terminal in neuron_tree.terminals:\n",
    "            if terminal.type==3:\n",
    "                continue\n",
    "            neuron_terminals[neuron]['terminals'].append(terminal.xyz)\n",
    "    np.save(data_path+'neuron_terminals.npy',neuron_terminals)\n",
    "else:\n",
    "    neuron_terminals = np.load(data_path+'neuron_terminals.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# which part to which amygdala"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brainstem_region_hierachy = {\n",
    "    'Cortex':['FRP','PL','ILA','ORB','ACA','AI','GU','VISC','PTLp','RSP','TEa','ECT','PERI','MOs','MOp','SSp','SSs','AUD','VIS'],\n",
    "    'OLF':['MOB','AOB','AON','TT','DP','PIR','NLOT','COA','PAA','TR'],\n",
    "    'HPF':['CA','DG','FC','IG','ENT','PAR','POST','PRE','SUB','ProS','HATA','APr'],\n",
    "    'Cortex-sp':['CLA','EP','LA','BLA','BMA','PA'],\n",
    "    'CNU':['STRd','ACB','FS','OT','LSX','AAA','CEA','MEA','PALd','PALv','PALm','PALc'],\n",
    "    'TH':['VAL','VM','VPL','VPM','VPLpc','VPMpc','PoT','SPA','SPFm','SPFp','PP','MG','LGd','LP','PO','POL','SGN','Eth','AV','AM','AD','IAM','IAD','LD',\n",
    "    'IMD','MD','SMT','PR','PVT','PT','RE','Xi','RH','CM','PCN','CL','PF','PIL','GENv','EPI','RT'],\n",
    "    'HY':['SO','ASO','PVH','PVa','PVi','ARH','ADP','AVP','AVPV','DMH','MEPO','MPO','OV','PD','PS','PVp','PVpo','SBPV','SCH','SFO','VMPO','VLPO','AHN',\n",
    "    'MBO','MPN','PMd','PMv','PVHd','VMH','PH','LHA','LPO','PST','PSTN','PeF','RCH','STN','TU','ZI','ME'],\n",
    "    'MB':['APN','NOT','NPC','PPT','RPF','AT','MPT','OP','VTA','SNr','SNc','SCs','SCm','NB','SAG','PBG','IC','RR','MRN','RN',\n",
    "    'PAG','PRC','INC','ND','Su3','CUN','PPN','IF','IPN','RL','CLI','DR','III','MA3','EW','VTN'],\n",
    "    'P':['PSV','PB','PCG','PG','PRNr','PRNc','SUT','TRN','V','P5','CS','LDT','NI','RPO','SOC','NLL','DTN','PDTg','SG','I5','SLC','SLD'],\n",
    "    'MY':['DCO','VCO','LAV','MV','SPIV','SUV','DMX','GR','CU','ECU','NTS','NTB','SPVC','SPVI','SPVO','Pa5','VII','GRN','IRN','LIN','LRNm','LRNp','MARN',\n",
    "    'MDRNd','MDRNv','PARN','PGRNd','PGRNl','PHY','VI','AMB','ICB','PAS','PRP','PPY','x','XII','y','IO','RM','RPA','RO'],\n",
    "    'CB':['VERM','HEM','FN','IP','DN','VeCB'],\n",
    "    'FT':['py','pyd'],\n",
    "    'SC':['unknow']\n",
    "}\n",
    "brain_mask = 0\n",
    "brain_big_region = ['Isocortex','OLF','HPF','CNU','HY','TH','MB','P','CB','MY']\n",
    "for i,j in enumerate(brain_big_region):\n",
    "    brain_mask = np.where(get_region_mask(j),i+1,brain_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_distribution = {}\n",
    "for i in neuron_information:\n",
    "    r = get_mask_value(brain_mask,np.array(neuron_information[i]['soma_position'])/10)-1\n",
    "    if r==-1:\n",
    "        continue\n",
    "    r = brain_big_region[r]\n",
    "    if r not in region_distribution:\n",
    "        region_distribution[r] = []\n",
    "    region_distribution[r].append(neuron_information[i]['axonal_projection_terminal_combined'])\n",
    "tmp = [[],[]]\n",
    "for i in brain_big_region:\n",
    "    if i not in region_distribution:\n",
    "        continue\n",
    "    a = np.sum(region_distribution[i],axis=0)\n",
    "    tmp[0].append(i)\n",
    "    tmp[1].append(a)\n",
    "#tmp[1] = np.array(tmp[1])/np.sum(tmp[1],axis=-1,keepdims=1)\n",
    "plt.matshow(tmp[1],cmap='hot')\n",
    "plt.xticks(range(len(amygdala_region)),amygdala_region)\n",
    "plt.yticks(range(len(tmp[0])),tmp[0])\n",
    "plt.title('Projection strength',fontsize=fontsize_title)\n",
    "plt.colorbar()\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(data_path+'projection_strength_to_AMY.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 0\n",
    "for i,j in zip(amygdala_region,np.array(tmp[1]).T):\n",
    "    k = j/np.max(j)\n",
    "    plt.bar(range(len(j)),k,bottom=c-0.5,label=i,width=0.5,color=color_pool[0])\n",
    "    plt.axhline(c-0.5,color='#AAAAAA')\n",
    "    c += 1\n",
    "plt.yticks(range(len(amygdala_region)),amygdala_region)\n",
    "plt.xticks(range(len(tmp[0])),tmp[0])\n",
    "plt.title('projection ratio')\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(data_path+'projection_ratio_to_AMY.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## projection density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain_mask = 0\n",
    "brain_big_region = ['FRP','PL','ILA','ORB','ACA','AI','GU','VISC','PTLp','RSP','TEa','ECT','PERI','MOs','MOp','SSp','SSs','AUD','VIS']\n",
    "for i,j in enumerate(brain_big_region):\n",
    "    brain_mask = np.where(get_region_mask(j),i+1,brain_mask)\n",
    "region_distribution = {}\n",
    "for i in neuron_information:\n",
    "    r = get_mask_value(brain_mask,np.array(neuron_information[i]['soma_position'])/10)-1\n",
    "    if r==-1:\n",
    "        continue\n",
    "    r = brain_big_region[r]\n",
    "    if r not in region_distribution:\n",
    "        region_distribution[r] = []\n",
    "    tmp = compute_flatten_position_and_depth([neuron_information[i]['soma_position']])[0]\n",
    "    if tmp[0]<=0:\n",
    "        continue\n",
    "    region_distribution[r].append({'name':i,'depth':tmp[0],'flatten':tmp[-1],'projection':neuron_information[i]['axonal_projection_terminal_combined']})\n",
    "tmp = {}\n",
    "for i in region_distribution:\n",
    "    for j in region_distribution[i]:\n",
    "        tmp[j['name']] = j\n",
    "region_distribution = tmp\n",
    "density_info = {}\n",
    "r = 250\n",
    "X = []\n",
    "for neuron in region_distribution:\n",
    "    if len(region_distribution[neuron]['flatten'])>0:\n",
    "        X.append(region_distribution[neuron]['flatten'])\n",
    "X = np.array(X)\n",
    "for neuron in region_distribution:\n",
    "    if len(region_distribution[neuron]['flatten'])>0:\n",
    "        tmp = np.linalg.norm(np.array(region_distribution[neuron]['flatten'])-X,axis=-1)\n",
    "        mask = np.array(tmp<r)\n",
    "        density_info[neuron] = np.sum(mask*(r-tmp))\n",
    "_max = np.max(list(density_info.values()))\n",
    "for neuron in region_distribution:\n",
    "    if len(region_distribution[neuron]['flatten'])>0:\n",
    "        region_distribution[neuron]['relative_density'] = density_info[neuron]/_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info = {}\n",
    "for index,t in enumerate(amygdala_region):\n",
    "    a = 0\n",
    "    for i in region_distribution:\n",
    "        a += 1\n",
    "        info[i] = {'soma_flatten_position':region_distribution[i]['flatten'],'projection_strength':np.sum(region_distribution[i]['projection'][index]),'relative_density':region_distribution[i]['relative_density'],}\n",
    "    plot_connection_density(info,save_path=data_path+f'projection_density_to_{t}_({a}_neurons).pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info = {}\n",
    "for index,t in enumerate(amygdala_region):\n",
    "    a = []\n",
    "    for i in region_distribution:\n",
    "        if np.sum(region_distribution[i]['projection'][index])>5:\n",
    "            a.append(i)\n",
    "    save_scene(a,color_pool[index],data_path+f'cortex_to_{t}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = {}\n",
    "for neuron in Scene.scene2List(data_path+'cortex_to_LA.nv'):\n",
    "    i = neuron['sampleid']+neuron['name']\n",
    "    r = get_mask_value(brain_mask,np.array(neuron_information[i]['soma_position'])/10)-1\n",
    "    if r==-1:\n",
    "        continue\n",
    "    r = brain_big_region[r]\n",
    "    if r not in A:\n",
    "        A[r] = []\n",
    "    A[r].append(i)\n",
    "for index,i in enumerate(A):\n",
    "    save_scene(A[i],color_pool[index],data_path+f'cortex_to_LA_{i}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Short and Long projection with AMY projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iso_mask = get_region_mask('Isocortex')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for neuron in neuron_information:\n",
    "    r = get_mask_value(brain_mask,np.array(neuron_information[i]['soma_position'])/10)-1\n",
    "    if r==-1:\n",
    "        continue\n",
    "    r = brain_big_region[r]\n",
    "    tmp = []\n",
    "    for terminal in neuron_terminals[neuron]['terminals']:\n",
    "        if if_inside_mask(Iso_mask,np.array(terminal)/10):\n",
    "            tmp.append(np.linalg.norm(neuron_information[i]['soma_position']-np.array(terminal)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for neuron in neuron_terminals:\n",
    "    depth = neuron_to_brainstem_information[neuron]['soma_depth']\n",
    "    if depth==0:\n",
    "        continue\n",
    "    tmp = []\n",
    "    for terminal in neuron_terminals[neuron]['terminals']:\n",
    "        if if_inside_mask(mask,np.array(terminal)/10):\n",
    "            tmp.append(terminal[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cortex->CEA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating isocortex data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CEA_data_path = data_path+'CEA/'\n",
    "create_dir(CEA_data_path)\n",
    "cortex_region = ['FRP','PL','ILA','ORB','ACA','AI','GU','VISC','PTLp','RSP','TEa','ECT','PERI','MOs','MOp','SSp','SSs','AUD','VIS']\n",
    "CEA_region = ['CEAl','CEAc','CEAm']\n",
    "cortex_region_to_subregion,cortex_subregion_to_region = get_subregions(cortex_region)\n",
    "if not os.path.exists(CEA_data_path+'Isocortex_neuron_info.npy'):\n",
    "    mask = get_region_mask('Isocortex')\n",
    "    Isocortex_neuron_info = {}\n",
    "    for i in amygdala_projection_neuron_name:\n",
    "        region = int(get_mask_value(mask,np.array(neuron_information[i]['soma_position'])/10))\n",
    "        if region>0:\n",
    "            Isocortex_neuron_info[i] = {}\n",
    "    for i in tool_bar(Isocortex_neuron_info):\n",
    "        neuron_property = iondata.getNeuronPropertyByID(i[:-7],i[-7:])\n",
    "        Isocortex_neuron_info[i]['soma_position'] = neuron_property['somapoint']\n",
    "        Isocortex_neuron_info[i]['soma_hemisphere'] = 'left' if neuron_property['somapoint'][2]<5700 else 'right'\n",
    "        a,b = region_layer_seperate(neuron_property['somaregion'])\n",
    "        Isocortex_neuron_info[i]['cortex_region'] = cortex_subregion_to_region[neuron_property['somaregion']]\n",
    "        Isocortex_neuron_info[i]['cortex_subregion'] = a[len(Isocortex_neuron_info[i]['cortex_region']):]\n",
    "        if Isocortex_neuron_info[i]['cortex_region']=='PTLp' and a[:4]!='PTLp':\n",
    "            Isocortex_neuron_info[i]['cortex_subregion'] = a[3:]\n",
    "        Isocortex_neuron_info[i]['cortex_layer'] = b\n",
    "        if len(Isocortex_neuron_info[i]['cortex_subregion'])!=0 and Isocortex_neuron_info[i]['cortex_subregion']=='-':\n",
    "            Isocortex_neuron_info[i]['cortex_subregion'] = Isocortex_neuron_info[i]['cortex_subregion'][1:]\n",
    "        Isocortex_neuron_info[i]['axonal_projection_length_combined'] = get_axonal_statics_by_regions(i,CEA_region,'projectregion')\n",
    "        Isocortex_neuron_info[i]['axonal_projection_length_left'] = get_axonal_statics_by_regions(i,CEA_region,'projectleftregion')\n",
    "        Isocortex_neuron_info[i]['axonal_projection_length_right'] = get_axonal_statics_by_regions(i,CEA_region,'projectrightregion')\n",
    "        Isocortex_neuron_info[i]['axonal_projection_terminal_combined'] = get_axonal_statics_by_regions(i,CEA_region,'terminalregion')\n",
    "        Isocortex_neuron_info[i]['axonal_projection_terminal_left'] = get_axonal_statics_by_regions(i,CEA_region,'terminalleftregion')\n",
    "        Isocortex_neuron_info[i]['axonal_projection_terminal_right'] = get_axonal_statics_by_regions(i,CEA_region,'terminalrightregion')\n",
    "        if Isocortex_neuron_info[i]['soma_hemisphere']=='right':\n",
    "            tmp = Isocortex_neuron_info[i]['axonal_projection_length_left'].copy()\n",
    "            Isocortex_neuron_info[i]['axonal_projection_length_left'] = Isocortex_neuron_info[i]['axonal_projection_length_right']\n",
    "            Isocortex_neuron_info[i]['axonal_projection_length_right'] = tmp.copy()\n",
    "            tmp = Isocortex_neuron_info[i]['axonal_projection_terminal_left'].copy()\n",
    "            Isocortex_neuron_info[i]['axonal_projection_terminal_left'] = Isocortex_neuron_info[i]['axonal_projection_terminal_right']\n",
    "            Isocortex_neuron_info[i]['axonal_projection_terminal_right'] = tmp.copy()\n",
    "        Isocortex_neuron_info[i]['dendrite_terminal_point'] = []\n",
    "        Isocortex_neuron_info[i]['axon_terminal_point'] = []\n",
    "        control_console_output(0)\n",
    "        neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "        control_console_output(1)\n",
    "        for terminal in neuron_tree.terminals:\n",
    "            if terminal.type==3:\n",
    "                Isocortex_neuron_info[i]['dendrite_terminal_point'].append(terminal.xyz)\n",
    "            else:\n",
    "                Isocortex_neuron_info[i]['axon_terminal_point'].append(terminal.xyz)\n",
    "        tmp = compute_flatten_position_and_depth([Isocortex_neuron_info[i]['soma_position']])[0]\n",
    "        Isocortex_neuron_info[i]['soma_depth'] = tmp[0]\n",
    "        Isocortex_neuron_info[i]['soma_outter_position'] = tmp[1]\n",
    "        Isocortex_neuron_info[i]['soma_inner_position'] = tmp[2]\n",
    "        Isocortex_neuron_info[i]['soma_flatten_position'] = tmp[3]\n",
    "    r = 50.0\n",
    "    X = []\n",
    "    for neuron in Isocortex_neuron_info:\n",
    "        X.append(Isocortex_neuron_info[neuron]['soma_inner_position'])\n",
    "    X = np.array(X)\n",
    "    for i in list(Isocortex_neuron_info.keys()):\n",
    "        if np.sum(Isocortex_neuron_info[i]['axonal_projection_terminal_combined'])==0:\n",
    "            del Isocortex_neuron_info[i]\n",
    "    density_info = {}\n",
    "    for neuron in Isocortex_neuron_info:\n",
    "        tmp = np.linalg.norm(np.array(Isocortex_neuron_info[neuron]['soma_inner_position'])-X,axis=-1)\n",
    "        mask = np.sum(tmp<r)\n",
    "        density_info[neuron] = np.sum(mask*(r-tmp))\n",
    "    _max = np.max(list(density_info.values()))\n",
    "    for neuron in Isocortex_neuron_info:\n",
    "        Isocortex_neuron_info[neuron]['relative_density'] = density_info[neuron]/_max\n",
    "    np.save(CEA_data_path+'Isocortex_neuron_info.npy',Isocortex_neuron_info)\n",
    "else:\n",
    "    Isocortex_neuron_info = np.load(CEA_data_path+'Isocortex_neuron_info.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Region statics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_distribution = {}\n",
    "for i in Isocortex_neuron_info:\n",
    "    r = Isocortex_neuron_info[i]['cortex_region']\n",
    "    if r not in region_distribution:\n",
    "        region_distribution[r] = []\n",
    "    region_distribution[r].append(i)\n",
    "for i,j in enumerate(list(region_distribution.keys())):\n",
    "    save_scene(region_distribution[j],color_pool[i%12],CEA_data_path+j)\n",
    "for i in region_distribution:\n",
    "    region_distribution[i] = len(region_distribution[i])\n",
    "tmp = list(region_distribution.keys())\n",
    "region_distribution['others'] = 0\n",
    "_sum = np.sum(list(region_distribution.values()))\n",
    "for i in tmp:\n",
    "    if (region_distribution[i]/_sum)<0.05:\n",
    "        region_distribution['others'] += region_distribution[i]\n",
    "        del region_distribution[i]\n",
    "plt.figure(figsize=(5,5))\n",
    "def func(pct, allvals):\n",
    "    absolute = int(np.round(pct/100.*np.sum(allvals)))\n",
    "    return f\"{pct:.1f}%\\n({absolute:d})\"\n",
    "plt.pie(region_distribution.values(),labels=region_distribution.keys(),radius=1,startangle=90,colors=color_pool[0:4],explode=[0.05 for i in region_distribution], autopct=lambda pct: func(pct,list(region_distribution.values())))\n",
    "plt.pie([10],colors=['w'],radius=0.4)\n",
    "plt.title('Region distribution',fontsize=fontsize_title)\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+'neurons_region_distribution.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Soma location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = []\n",
    "Y = []\n",
    "for i in Isocortex_neuron_info:\n",
    "    if np.sum(Isocortex_neuron_info[i]['axonal_projection_terminal_combined'])>0:\n",
    "        X.append(Isocortex_neuron_info[i]['soma_position'])\n",
    "        if X[-1][2]>5700:\n",
    "            X[-1][2] = 11400-X[-1][2]\n",
    "        Y.append([1,0,0])\n",
    "neuronvis = nv.neuronVis(size=(1920,1080))\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "draw_points(neuronvis,X,Y)\n",
    "neuronvis.render.setPointSize(8)\n",
    "set_view(neuronvis,'left')\n",
    "neuronvis.render.savepng(CEA_data_path+'soma_location.png')\n",
    "neuronvis.render.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connection density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info = {}\n",
    "for i in Isocortex_neuron_info:\n",
    "    r = Isocortex_neuron_info[i]['cortex_region']\n",
    "    info[i] = {'soma_flatten_position':Isocortex_neuron_info[i]['soma_flatten_position'],'projection_strength':np.sum(Isocortex_neuron_info[i]['axonal_projection_terminal_combined']),'relative_density':Isocortex_neuron_info[i]['relative_density'],}\n",
    "plot_connection_density(info,save_path=CEA_data_path+'connection_density.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for r in CEA_region:\n",
    "    info = {}\n",
    "    for i in Isocortex_neuron_info:\n",
    "        info[i] = {'soma_flatten_position':Isocortex_neuron_info[i]['soma_flatten_position'],'projection_strength':Isocortex_neuron_info[i]['axonal_projection_terminal_combined'][CEA_region.index(r)],'relative_density':Isocortex_neuron_info[i]['relative_density'],}\n",
    "    plot_connection_density(info,save_path=CEA_data_path+'connection_density_'+r+'.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coprojection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chi2_contingency\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "_flatmap = cv.imread('../../resource/flatmapedge2.tif',1)\n",
    "def cluster_neuron(neurons,regions,X,save_path,circle_size=20,color_threshold=5):\n",
    "    #X_length = np.linalg.norm(X,axis=-1)\n",
    "    #X_length = np.max(X,axis=-1)\n",
    "    #X_length = np.sum(X,axis=-1)\n",
    "    #feature = X/np.where(X_length==0,1,X_length).reshape(-1,1)+1e-8\n",
    "    #a = pdist(feature,metric='cosine')\n",
    "    feature = np.array(X)*1\n",
    "    print(feature.shape)\n",
    "    feature[:,:2] = feature[:,:2]**0.5\n",
    "    feature[:,feature.shape[1]//2:feature.shape[1]//2+2] = feature[:,feature.shape[1]//2:feature.shape[1]//2+2]**0.5\n",
    "    feature_length = np.max(feature,axis=-1)\n",
    "    feature = feature/np.where(feature_length==0,1,feature_length).reshape(-1,1)+1e-8\n",
    "    a = feature\n",
    "    fig = plt.figure(figsize=(7,5.5))\n",
    "    ax_dendro = fig.add_axes([0.1,0.8,0.39,0.2])\n",
    "    cortex_cluster_object = hierarchy.dendrogram(hierarchy.linkage(a,'ward'),distance_sort=True,orientation='top',ax=ax_dendro,color_threshold=color_threshold)\n",
    "    ax_dendro.set_xticks([],[])\n",
    "    ax_dendro.set_yticks([],[])\n",
    "    ax_dendro.set_title('Clustering',fontsize=fontsize_label)\n",
    "    cluster_object_parse = [[],[]]\n",
    "    neuron_to_cluster = {}\n",
    "    for i,j in zip(cortex_cluster_object['leaves'],cortex_cluster_object['leaves_color_list']):\n",
    "        neuron = neurons[i]\n",
    "        projection_matrix = feature[i]\n",
    "        try:\n",
    "            cluster_object_parse[0].append(color_pool_hex.index(j))\n",
    "            cluster_object_parse[1].append(projection_matrix)\n",
    "            neuron_to_cluster[neuron] = color_pool_hex.index(j)\n",
    "        except:\n",
    "            cluster_object_parse[0].append(-1)\n",
    "            cluster_object_parse[1].append(projection_matrix)\n",
    "            neuron_to_cluster[neuron] = -1\n",
    "    cluster_object_parse[1] = np.array(cluster_object_parse[1])\n",
    "    ax_matrix = fig.add_axes([0.1,0.0,0.39,0.8])\n",
    "    ax_matrix.matshow(cluster_object_parse[1].T,cmap='BuGn',aspect='auto')\n",
    "    _ = ax_matrix.set_xticks([],[])\n",
    "    ax_matrix.set_xlabel('neurons = '+str(len(neurons)),fontsize=fontsize_label)\n",
    "    cluster_object_parse_visualization_yticks = [[],[]]\n",
    "    cluster_object_parse_visualization_yticks[0] = [i for i in range(len(regions))]\n",
    "    cluster_object_parse_visualization_yticks[1] = [i if i!='unknow' else 'Spine' for i in regions]\n",
    "    #_ = ax_matrix.set_yticks(cluster_object_parse_visualization_yticks[0],cluster_object_parse_visualization_yticks[1],fontsize=fontsize_ticks)\n",
    "    ax_matrix.set_yticks(cluster_object_parse_visualization_yticks[0])\n",
    "    ax_matrix.set_yticklabels(cluster_object_parse_visualization_yticks[1])\n",
    "    ax_matrix.set_ylabel('Traget',fontsize=fontsize_label,rotation=-90).set_verticalalignment('top')\n",
    "    ax_soma = fig.add_axes([0.5,0.5,0.4,0.5])\n",
    "    cluster_info = {}\n",
    "    for i,j in enumerate(neurons):\n",
    "        _cluster = neuron_to_cluster[j]\n",
    "        if _cluster==-1:\n",
    "            print(j)\n",
    "            continue\n",
    "        if _cluster not in cluster_info:\n",
    "            cluster_info[_cluster] = []\n",
    "        cluster_info[_cluster].append([j,Isocortex_neuron_info[j]['soma_flatten_position']])\n",
    "    flatmap_tmp = _flatmap*1\n",
    "    for i in cluster_info.keys():\n",
    "        tmp3 = []\n",
    "        for j in cluster_info[i]:\n",
    "            if np.sum(abs(np.array(j[1])))!=0:\n",
    "                cv.circle(flatmap_tmp,[int(j[1][0])*2, int(j[1][1])*2],circle_size,color_pool[i]*255,-1)\n",
    "            tmp3.append(j[0])\n",
    "        save_scene(tmp3,color_pool[i],save_path+'_cluster_'+str(i+1))\n",
    "    ax_soma.imshow(flatmap_tmp)\n",
    "    ax_soma.imshow(flatmap,zorder=1)\n",
    "    ax_soma.set_xlim(flatten_map_r[0],flatten_map_r[1])\n",
    "    ax_soma.set_ylim(flatten_map_r[2],flatten_map_r[3])\n",
    "    ax_soma.set_axis_off()\n",
    "    ax_soma.set_title('Soma distribution',fontsize=fontsize_label)\n",
    "    cluster_info_transform = {}\n",
    "    for i in cluster_info:\n",
    "        for j in cluster_info[i]:\n",
    "            cluster_info_transform[j[0]] = i\n",
    "    ax_feature = fig.add_axes([0.5,0.0,0.4,0.45])\n",
    "    r = {}\n",
    "    Sp = {}\n",
    "    eee = {}\n",
    "    for i in cluster_info.keys():\n",
    "        eee[i] = []\n",
    "        r[i] = []\n",
    "        for j in cluster_info[i]:\n",
    "            r[i].append(X[neurons.index(j[0])])\n",
    "        tmp = np.mean(r[i],axis=0)\n",
    "        eee[i] = np.mean(eee[i])/np.where(np.max(tmp)==0,1,np.max(tmp))\n",
    "        tmp /= np.where(np.sum(tmp)==0,1,np.sum(tmp))\n",
    "        Sp[i] = tmp\n",
    "    _cls_list = sorted(list(cluster_info.keys()))\n",
    "    ck = [np.argmax(Sp[i]) for i in _cls_list]\n",
    "    for index,i in enumerate(sorted(range(len(ck)),key=lambda k:ck[k])):\n",
    "        ax_feature.barh([-j for j in range(len(Sp[_cls_list[i]]))],Sp[_cls_list[i]].tolist(),left=index,color=color_pool[_cls_list[i]])\n",
    "    zz = np.zeros((len(_cls_list),len(_cls_list)))\n",
    "    for i in range(len(_cls_list)):\n",
    "        for j in range(-1,i):\n",
    "            op = np.array([Sp[_cls_list[i]],Sp[_cls_list[j]]])\n",
    "            op += np.where(np.sum(op,axis=0).reshape(1,-1)==0,1,np.sum(op,axis=0).reshape(1,-1)) # zeor expectation case\n",
    "            zz[i,j] = chi2_contingency(op)[1]\n",
    "            zz[j,i] = -1\n",
    "    for i in range(len(cluster_info)-1):\n",
    "        ax_feature.axvline(i+1,linewidth=0.5,color='#AAAAAA')\n",
    "    ax_feature.set_xlim(0,len(cluster_info))\n",
    "    ax_feature.set_xticks([],[])\n",
    "    #ax_feature.set_yticks([-j for j in range(len(tmp))],cluster_object_parse_visualization_yticks[1],fontsize=fontsize_ticks if len(cluster_info)<13 else fontsize_ticks*0.4)\n",
    "    ax_feature.set_yticks([-j for j in range(len(tmp))])\n",
    "    ax_feature.set_yticklabels(cluster_object_parse_visualization_yticks[1])\n",
    "    ax_feature.yaxis.tick_right()\n",
    "    ax_feature.set_xlabel('Clusters',fontsize=fontsize_label)\n",
    "    ax_feature.set_title('Projection feature',fontsize=fontsize_label)\n",
    "    plt.rcParams['pdf.fonttype'] = 42\n",
    "    plt.rcParams['ps.fonttype'] = 42\n",
    "    plt.savefig(save_path+'_overall.pdf',format='pdf',bbox_inches='tight',dpi=100)\n",
    "    plt.show()\n",
    "    return cluster_info,zz,Sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Coprojection_region = ['CTX','CNU','HY','TH','MB','P','CB','MY']\n",
    "if not os.path.exists(CEA_data_path+'Coprojection_info.npy'):\n",
    "    Coprojection_info = {}\n",
    "    for i in Isocortex_neuron_info:\n",
    "        Coprojection_info[i] = {}\n",
    "        Coprojection_info[i]['axonal_projection_length_combined'] = get_axonal_statics_by_regions(i,Coprojection_region,'projectregion')\n",
    "        Coprojection_info[i]['axonal_projection_length_left'] = get_axonal_statics_by_regions(i,Coprojection_region,'projectleftregion')\n",
    "        Coprojection_info[i]['axonal_projection_length_right'] = get_axonal_statics_by_regions(i,Coprojection_region,'projectrightregion')\n",
    "        Coprojection_info[i]['axonal_projection_terminal_combined'] = get_axonal_statics_by_regions(i,Coprojection_region,'terminalregion')\n",
    "        Coprojection_info[i]['axonal_projection_terminal_left'] = get_axonal_statics_by_regions(i,Coprojection_region,'terminalleftregion')\n",
    "        Coprojection_info[i]['axonal_projection_terminal_right'] = get_axonal_statics_by_regions(i,Coprojection_region,'terminalrightregion')\n",
    "        if Isocortex_neuron_info[i]['soma_hemisphere']=='right':\n",
    "            tmp = Coprojection_info[i]['axonal_projection_length_left'].copy()\n",
    "            Coprojection_info[i]['axonal_projection_length_left'] = Coprojection_info[i]['axonal_projection_length_right']\n",
    "            Coprojection_info[i]['axonal_projection_length_right'] = tmp.copy()\n",
    "            tmp = Coprojection_info[i]['axonal_projection_terminal_left'].copy()\n",
    "            Coprojection_info[i]['axonal_projection_terminal_left'] = Coprojection_info[i]['axonal_projection_terminal_right']\n",
    "            Coprojection_info[i]['axonal_projection_terminal_right'] = tmp.copy()\n",
    "    np.save(CEA_data_path+'Coprojection_info.npy',Coprojection_info)\n",
    "else:\n",
    "    Coprojection_info = np.load(CEA_data_path+'Coprojection_info.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brainstem_region_hierachy = {\n",
    "    'Cortex':['FRP','PL','ILA','ORB','ACA','AI','GU','VISC','PTLp','RSP','TEa','ECT','PERI','MOs','MOp','SSp','SSs','AUD','VIS'],\n",
    "    'OLF':['MOB','AOB','AON','TT','DP','PIR','NLOT','COA','PAA','TR'],\n",
    "    'HPF':['CA','DG','FC','IG','ENT','PAR','POST','PRE','SUB','ProS','HATA','APr'],\n",
    "    'Cortex-sp':['CLA','EP','LA','BLA','BMA','PA'],\n",
    "    'CNU':['STRd','ACB','FS','OT','LSX','AAA','CEA','MEA','PALd','PALv','PALm','PALc'],\n",
    "    'TH':['VAL','VM','VPL','VPM','VPLpc','VPMpc','PoT','SPA','SPFm','SPFp','PP','MG','LGd','LP','PO','POL','SGN','Eth','AV','AM','AD','IAM','IAD','LD',\n",
    "    'IMD','MD','SMT','PR','PVT','PT','RE','Xi','RH','CM','PCN','CL','PF','PIL','GENv','EPI','RT'],\n",
    "    'HY':['SO','ASO','PVH','PVa','PVi','ARH','ADP','AVP','AVPV','DMH','MEPO','MPO','OV','PD','PS','PVp','PVpo','SBPV','SCH','SFO','VMPO','VLPO','AHN',\n",
    "    'MBO','MPN','PMd','PMv','PVHd','VMH','PH','LHA','LPO','PST','PSTN','PeF','RCH','STN','TU','ZI','ME'],\n",
    "    'MB':['APN','NOT','NPC','PPT','RPF','AT','MPT','OP','VTA','SNr','SNc','SCs','SCm','NB','SAG','PBG','IC','RR','MRN','RN',\n",
    "    'PAG','PRC','INC','ND','Su3','CUN','PPN','IF','IPN','RL','CLI','DR','III','MA3','EW','VTN'],\n",
    "    'P':['PSV','PB','PCG','PG','PRNr','PRNc','SUT','TRN','V','P5','CS','LDT','NI','RPO','SOC','NLL','DTN','PDTg','SG','I5','SLC','SLD'],\n",
    "    'MY':['DCO','VCO','LAV','MV','SPIV','SUV','DMX','GR','CU','ECU','NTS','NTB','SPVC','SPVI','SPVO','Pa5','VII','GRN','IRN','LIN','LRNm','LRNp','MARN',\n",
    "    'MDRNd','MDRNv','PARN','PGRNd','PGRNl','PHY','VI','AMB','ICB','PAS','PRP','PPY','x','XII','y','IO','RM','RPA','RO'],\n",
    "    'CB':['VERM','HEM','FN','IP','DN','VeCB'],\n",
    "    'FT':['py','pyd'],\n",
    "    'SC':['unknow']\n",
    "}\n",
    "cortex_region_to_subregion,cortex_subregion_to_region = get_subregions(cortex_region)\n",
    "brainstem_region = []\n",
    "for a in brainstem_region_hierachy:\n",
    "    for b in brainstem_region_hierachy[a]:\n",
    "        brainstem_region.append(b)\n",
    "brainstem_region_to_subregion,brainstem_subregion_to_region = get_subregions(brainstem_region)\n",
    "if not os.path.exists(CEA_data_path+'Coprojection_detailed_info.npy'):\n",
    "    Coprojection_detailed_info = {}\n",
    "    for i in Isocortex_neuron_info:\n",
    "        Coprojection_detailed_info[i] = {}\n",
    "        Coprojection_detailed_info[i]['axonal_projection_length_combined'] = get_axonal_statics_by_regions(i,brainstem_region,'projectregion')\n",
    "        Coprojection_detailed_info[i]['axonal_projection_length_left'] = get_axonal_statics_by_regions(i,brainstem_region,'projectleftregion')\n",
    "        Coprojection_detailed_info[i]['axonal_projection_length_right'] = get_axonal_statics_by_regions(i,brainstem_region,'projectrightregion')\n",
    "        Coprojection_detailed_info[i]['axonal_projection_terminal_combined'] = get_axonal_statics_by_regions(i,brainstem_region,'terminalregion')\n",
    "        Coprojection_detailed_info[i]['axonal_projection_terminal_left'] = get_axonal_statics_by_regions(i,brainstem_region,'terminalleftregion')\n",
    "        Coprojection_detailed_info[i]['axonal_projection_terminal_right'] = get_axonal_statics_by_regions(i,brainstem_region,'terminalrightregion')\n",
    "        if Isocortex_neuron_info[i]['soma_hemisphere']=='right':\n",
    "            tmp = Coprojection_detailed_info[i]['axonal_projection_length_left'].copy()\n",
    "            Coprojection_detailed_info[i]['axonal_projection_length_left'] = Coprojection_detailed_info[i]['axonal_projection_length_right']\n",
    "            Coprojection_detailed_info[i]['axonal_projection_length_right'] = tmp.copy()\n",
    "            tmp = Coprojection_detailed_info[i]['axonal_projection_terminal_left'].copy()\n",
    "            Coprojection_detailed_info[i]['axonal_projection_terminal_left'] = Coprojection_detailed_info[i]['axonal_projection_terminal_right']\n",
    "            Coprojection_detailed_info[i]['axonal_projection_terminal_right'] = tmp.copy()\n",
    "    np.save(CEA_data_path+'Coprojection_detailed_info.npy',Coprojection_detailed_info)\n",
    "else:\n",
    "    Coprojection_detailed_info = np.load(CEA_data_path+'Coprojection_detailed_info.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 0.05\n",
    "X = []\n",
    "neuron_list = list(Isocortex_neuron_info.keys())\n",
    "target = Coprojection_region+[i+'_c' for i in Coprojection_region]\n",
    "for i in neuron_list:\n",
    "    X.append(Coprojection_info[i]['axonal_projection_terminal_left']+Coprojection_info[i]['axonal_projection_terminal_right'])\n",
    "E,zz,debug = cluster_neuron(neuron_list,target,X,save_path=CEA_data_path,circle_size=15,color_threshold=5.5)\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "plt.matshow(zz>0.05,vmin=0,vmax=1,cmap='GnBu')\n",
    "plt.title('Dependency',fontsize=fontsize_title)\n",
    "plt.savefig(CEA_data_path+'cluster_dependency.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()\n",
    "for i in debug:\n",
    "    _max = np.max(debug[i])\n",
    "    c = []\n",
    "    for j,k in enumerate(debug[i]):\n",
    "        if k>_max*0.5:\n",
    "            c.append(target[j])\n",
    "    print(i,c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### coprojection analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = {}\n",
    "for cluster in [8,7,9,10,5,1,6,4,3,2]:\n",
    "    X[cluster] = []\n",
    "    for i in Scene.scene2List(CEA_data_path+f'_cluster_{cluster}.nv'):\n",
    "        neuron = i['sampleid']+i['name']\n",
    "        X[cluster].append(Isocortex_neuron_info[neuron]['axonal_projection_terminal_left']+Isocortex_neuron_info[neuron]['axonal_projection_terminal_right'])\n",
    "fig, ax = plt.subplots(1,len(X),figsize=(2*len(X),2))\n",
    "for i,j in enumerate(X):\n",
    "    tmp = np.mean(np.array(X[j]),axis=0)\n",
    "    ax[i].barh(range(3,0,-1),tmp[:3],color=color_pool[0],height=0.4)\n",
    "    ax[i].barh(range(3,0,-1),tmp[3:],left=tmp[:3],color=color_pool[3],height=0.4)\n",
    "    ax[i].set_title(f'cluster {j}')\n",
    "    if i==0:\n",
    "        ax[i].set_yticks(range(3,0,-1),CEA_region)\n",
    "    else:\n",
    "        ax[i].set_axis_off()\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+'cluster_in_CEA_projection.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### coprojection differences in CEA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_color = {8:0,7:0,9:3,10:3,5:1,1:1,6:1,4:1,3:1,2:1}\n",
    "X = []\n",
    "Y = []\n",
    "#mask = get_region_mask('CEA')\n",
    "for cluster in _color:\n",
    "    for i in Scene.scene2List(CEA_data_path+f'_cluster_{cluster}.nv'):\n",
    "        neuron = i['sampleid']+i['name']\n",
    "        for terminal in Isocortex_neuron_info[neuron]['axon_terminal_point']:\n",
    "            if if_inside_mask(mask,np.array(terminal)/10):\n",
    "                X.append(terminal.copy())\n",
    "                if X[-1][2]>5700:\n",
    "                    X[-1][2] = 11400-X[-1][2]\n",
    "                Y.append(color_pool[_color[cluster]])\n",
    "X = np.array(X)\n",
    "Y = np.array(Y)\n",
    "neuronvis = nv.neuronVis(size=(1920,1080))\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "draw_points(neuronvis,X,Y.tolist())\n",
    "neuronvis.render.setPointSize(2)\n",
    "neuronvis.addRegion('CEA')\n",
    "set_view(neuronvis,'dorsal')\n",
    "neuronvis.render.savepng(CEA_data_path+'cluster_in_CEA_projection_3D.png')\n",
    "neuronvis.render.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,figsize=(8,5))\n",
    "ax[0].scatter(X[:,2],-X[:,0],color=Y,s=1)\n",
    "ax[0].set_axis_off()\n",
    "ax[1].scatter(X[:,0],-X[:,1],color=Y,s=1)\n",
    "ax[1].set_axis_off()\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+'cluster_in_CEA_projection_2D_xyz.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in _color:\n",
    "    mask = np.sum((Y-color_pool[i-1])**2,axis=-1)==0\n",
    "    fig, ax = plt.subplots(1,2,figsize=(8,5))\n",
    "    ax[0].scatter(X[:,2][mask],-X[:,0][mask],color=Y[mask],s=1)\n",
    "    ax[0].set_axis_off()\n",
    "    ax[1].scatter(X[:,0][mask],-X[:,1][mask],color=Y[mask],s=1)\n",
    "    ax[1].set_axis_off()\n",
    "    plt.rcParams['pdf.fonttype'] = 42\n",
    "    plt.rcParams['ps.fonttype'] = 42\n",
    "    plt.savefig(CEA_data_path+f'cluster_in_CEA_projection_2D_xyz_cluster_{i}.pdf',format='pdf',bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = []\n",
    "Y = []\n",
    "mask = get_region_mask('CEA')\n",
    "for cluster in [8,7]:\n",
    "    for i in Scene.scene2List(CEA_data_path+f'_cluster_{cluster}.nv'):\n",
    "        neuron = i['sampleid']+i['name']\n",
    "        for terminal in Isocortex_neuron_info[neuron]['axon_terminal_point']:\n",
    "            if if_inside_mask(mask,np.array(terminal)/10):\n",
    "                X.append(terminal.copy())\n",
    "                if X[-1][2]>5700:\n",
    "                    X[-1][2] = 11400-X[-1][2]\n",
    "                Y.append(color_pool[cluster-7])\n",
    "X = np.array(X)\n",
    "fig, ax = plt.subplots(1,2,figsize=(8,5))\n",
    "ax[0].scatter(X[:,2],-X[:,0],color=Y,s=1)\n",
    "ax[0].set_axis_off()\n",
    "ax[1].scatter(X[:,0],-X[:,1],color=Y,s=1)\n",
    "ax[1].set_axis_off()\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+f'cluster_in_CEA_projection_2D_xyz_cluster_8_7.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### each cluster represents something aspect of CEA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_info = {}\n",
    "for ai,cluster in enumerate([8,7,9,10,5,1,6,4,3,2]):\n",
    "    cluster_info[cluster] = {'projection':[],'weights':[]}\n",
    "    for i in Scene.scene2List(CEA_data_path+f'_cluster_{cluster}.nv'):\n",
    "        neuron = i['sampleid']+i['name']\n",
    "        cluster_info[cluster]['projection'].append(Coprojection_detailed_info[neuron]['axonal_projection_terminal_combined'])\n",
    "        cluster_info[cluster]['weights'].append(Isocortex_neuron_info[neuron]['axonal_projection_terminal_combined'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "tmp = 0\n",
    "for ai,cluster in enumerate([8,7,9,10,5,1,6,4,3,2]):\n",
    "    a = np.mean(np.array(cluster_info[cluster]['projection'])*(np.sum(cluster_info[cluster]['weights'],axis=-1).reshape(-1,1)),axis=0)\n",
    "    a /= np.max(a)\n",
    "    plt.bar(range(len(a)),a,bottom=9-ai,color=color_pool[_color[cluster]])\n",
    "    tmp = np.where(a>tmp,a,tmp)\n",
    "xticks = [[],[]]\n",
    "for i in range(len(tmp)):\n",
    "    if tmp[i]>0.3:\n",
    "        xticks[0].append(i)\n",
    "        xticks[1].append(brainstem_region[i])\n",
    "plt.xticks(xticks[0],xticks[1],rotation=-90)\n",
    "plt.yticks([i+0.5 for i in range(10)],[8,7,9,10,5,1,6,4,3,2][::-1])\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+'cluster_detailed_projection.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neuronvis = nv.neuronVis(size=(1200,900))\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "neuronvis.render.setPointSize(1)\n",
    "for cluster in E:\n",
    "    neuronvis.clear(root=False)\n",
    "    X = []\n",
    "    for neuron in [i[0] for i in E[cluster]]:\n",
    "        if neuron in ['195617014.swc']:\n",
    "            continue\n",
    "        neuron_tree = iondata.getNeuronTreeByID(neuron[:-7],neuron[-7:])\n",
    "        neuron_tree.somaRadius = 10\n",
    "        neuronvis.addNeuronTree(neuron_tree,axonColor=color_pool[cluster],somaColor=[0.5,0.5,0.5],dendriteColor=[0.5,0.5,0.5],mirrorToRight=True,isLine=True,dendriteHide=True)\n",
    "        for terminal in neuron_tree.terminals:\n",
    "            X.append(terminal.xyz.copy())\n",
    "            if neuron_tree.root.z<7500:\n",
    "                X[-1][2] = 11400-X[-1][2]\n",
    "    for i in ['anterior','posterior','left','right','dorsal','ventral']:\n",
    "        set_view(neuronvis,i)\n",
    "        neuronvis.render.savepng(CEA_data_path+f'cluster_{cluster+1}_axon_{i}_.png')\n",
    "    neuronvis.clear(root=False)\n",
    "    draw_points(neuronvis,X,[color_pool[cluster].tolist() for i in X])\n",
    "    for i in ['anterior','posterior','left','right','dorsal','ventral']:\n",
    "        set_view(neuronvis,i)\n",
    "        neuronvis.render.savepng(CEA_data_path+f'cluster_{cluster+1}_terminal_{i}_.png')\n",
    "neuronvis.render.closeWindow()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Laminar pattern"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### data generating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 5\n",
    "lamina_info_CEA = {}\n",
    "lamina_info_non_CEA = {}\n",
    "for i in Isocortex_neuron_info:\n",
    "    region = Isocortex_neuron_info[i]['cortex_region']\n",
    "    if region not in ['AI']:\n",
    "        continue\n",
    "    layer = Isocortex_neuron_info[i]['cortex_layer']\n",
    "    if region=='CLA':\n",
    "        print(1)\n",
    "        layer = 7\n",
    "    if np.sum(Isocortex_neuron_info[i]['axonal_projection_terminal_combined'])<threshold:\n",
    "        g = lamina_info_non_CEA\n",
    "    else:\n",
    "        g = lamina_info_CEA\n",
    "    if layer not in g:\n",
    "        g[layer] = []\n",
    "    g[layer].append(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### cell number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5,3))\n",
    "x_ticks = ['1','2/3','5','6']\n",
    "y = []\n",
    "for i in x_ticks:\n",
    "    y.append(len(lamina_info_CEA[i]))\n",
    "plt.bar(range(len(x_ticks)),y,width=0.75)\n",
    "plt.xticks(range(len(x_ticks)),x_ticks)\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+'lamina_cell_number.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5,3))\n",
    "x_ticks = ['1','2/3','5','6']\n",
    "y = []\n",
    "for i in x_ticks:\n",
    "    y.append(len(lamina_info_CEA[i])/(len(lamina_info_CEA[i])+len(lamina_info_non_CEA[i])))\n",
    "plt.bar(range(len(x_ticks)),y,width=0.75)\n",
    "plt.xticks(range(len(x_ticks)),x_ticks)\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+'lamina_cell_ratio.pdf',format='pdf',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### coprojection in cortex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamina = ['2/3','5','6']\n",
    "lamina_color = {'2/3':0,'5':3,'6':1}\n",
    "mask = get_region_mask('Isocortex')\n",
    "X = {}\n",
    "Y = {}\n",
    "for i in lamina:\n",
    "    X[i] = []\n",
    "    Y[i] = []\n",
    "    for j in lamina_info_CEA[i]:\n",
    "        for k in Isocortex_neuron_info[j]['axon_terminal_point']:\n",
    "            if get_mask_value(mask,np.array(k)/10):\n",
    "                X[i].append(k.copy())\n",
    "                Y[i].append(color_pool[lamina_color[i]].tolist())\n",
    "X2 = {}\n",
    "for i in X:\n",
    "    X2[i] = {'depth':[],'flatten':[]}\n",
    "    tmp = compute_flatten_position_and_depth(X[i])\n",
    "    tmp2 = [[],[]]\n",
    "    for j in tmp:\n",
    "        X2[i]['depth'].append(j[0])\n",
    "        X2[i]['flatten'].append(j[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "for i,j in enumerate(lamina):\n",
    "    ax = plt.subplot(1,len(lamina),i+1)\n",
    "    tmp = flatmap2[...,:3]*1\n",
    "    color = color_pool[lamina_color[j]]#*255\n",
    "    color = color.astype(float)\n",
    "    for k in X2[j]['flatten']:\n",
    "        cv.circle(tmp,[int(np.round(k[0]))*2, int(np.round(k[1]))*2],5,color,-1)\n",
    "    ax.matshow(tmp)\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_xticks([],[])\n",
    "    ax.set_yticks([],[])\n",
    "plt.suptitle('Coprojection in cortex')\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+f'Coprojection_in_cortex_terminal_lamina.svg',format='svg',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "for i,j in enumerate(lamina):\n",
    "    ax = plt.subplot(1,len(lamina),i+1)\n",
    "    tmp = flatmap2[...,:3]*1\n",
    "    color = color_pool[lamina_color[j]]#*255\n",
    "    color = color.astype(float)\n",
    "    a = []\n",
    "    for k in X2[j]['depth']:\n",
    "        if k==0:\n",
    "            continue\n",
    "        a.append(k)\n",
    "    ax.hist(a,bins=200,range=[0,30000],density=True)\n",
    "plt.suptitle('Coprojection in cortex')\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+f'Coprojection_in_cortex_terminal_lamina_depth',format='pdf',bbox_inches='tight',dpi=500)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "for i,j in enumerate(lamina):\n",
    "    ax = plt.subplot(1,len(lamina),i+1)\n",
    "    tmp = flatmap2[...,:3]*1\n",
    "    color = color_pool[lamina_color[j]]#*255\n",
    "    color = color.astype(float)\n",
    "    a = []\n",
    "    for k,l in zip(X2[j]['flatten'],X2[j]['depth']):\n",
    "        a.append([k[0],l])\n",
    "    a = np.array(a)\n",
    "    ax.scatter(a[:,0],-a[:,1],s=0.5,color=color)\n",
    "plt.suptitle('Coprojection in cortex')\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+f'Coprojection_in_cortex_terminal_lamina_depth_visualization',format='pdf',bbox_inches='tight',dpi=500)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neuronvis = nv.neuronVis(size=(1200,900))\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "neuronvis.render.setPointSize(1)\n",
    "for cluster in E:\n",
    "    neuronvis.clear(root=False)\n",
    "    X = []\n",
    "    for neuron in [i[0] for i in E[cluster]]:\n",
    "        if neuron in ['195617014.swc']:\n",
    "            continue\n",
    "        neuron_tree = iondata.getNeuronTreeByID(neuron[:-7],neuron[-7:])\n",
    "        neuron_tree.somaRadius = 10\n",
    "        neuronvis.addNeuronTree(neuron_tree,axonColor=color_pool[cluster],somaColor=[0.5,0.5,0.5],dendriteColor=[0.5,0.5,0.5],mirrorToRight=True,isLine=True,dendriteHide=True)\n",
    "        for terminal in neuron_tree.terminals:\n",
    "            X.append(terminal.xyz.copy())\n",
    "            if neuron_tree.root.z<7500:\n",
    "                X[-1][2] = 11400-X[-1][2]\n",
    "    for i in ['anterior','posterior','left','right','dorsal','ventral']:\n",
    "        set_view(neuronvis,i)\n",
    "        neuronvis.render.savepng(CEA_data_path+f'cluster_{cluster+1}_axon_{i}_.png')\n",
    "    neuronvis.clear(root=False)\n",
    "    draw_points(neuronvis,X,[color_pool[cluster].tolist() for i in X])\n",
    "    for i in ['anterior','posterior','left','right','dorsal','ventral']:\n",
    "        set_view(neuronvis,i)\n",
    "        neuronvis.render.savepng(CEA_data_path+f'cluster_{cluster+1}_terminal_{i}_.png')\n",
    "neuronvis.render.closeWindow()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### terminal in CEA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = get_region_mask('CEA')\n",
    "X = {}\n",
    "Y = {}\n",
    "for i in lamina_info_CEA:\n",
    "    if i not in lamina_color:\n",
    "        continue\n",
    "    X[i] = []\n",
    "    Y[i] = []\n",
    "    for j in lamina_info_CEA[i]:\n",
    "        for k in Isocortex_neuron_info[j]['axon_terminal_point']:\n",
    "            if get_mask_value(mask,np.array(k)/10):\n",
    "                X[i].append(k.copy())\n",
    "                Y[i].append(color_pool[lamina_color[i]].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "for i,j in enumerate(lamina_color):\n",
    "    tmp = np.array(X[j])\n",
    "    ax = plt.subplot(1,len(lamina),i+1)\n",
    "    ax.scatter(tmp[:,2],tmp[:,1],color=Y[j],s=0.01)\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_xticks([],[])\n",
    "    ax.set_yticks([],[])\n",
    "plt.suptitle('Coprojection in cortex')\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "plt.savefig(CEA_data_path+f'Coprojection_in_CEA_terminal_lamina.svg',format='svg',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
